function [ys,params,check] = Ramsey_1s2c_perfect2_steadystate(ys,exo,M_,options_)
% function [ys,params,check] = rbc_optimal_tax_soe_Ramsey2_steadystate(ys,exo,M_,options_)
% computes the steady state for the rbc_optimal_tax_soe_Ramsey2.mod and uses a numerical
% solver to do so
% Inputs: 
%   - ys        [vector] vector of initial values for the steady state of
%                   the endogenous variables
%   - exo       [vector] vector of values for the exogenous variables
%   - M_        [structure] Dynare model structure
%   - options   [structure] Dynare options structure
%
% Output: 
%   - ys        [vector] vector of steady state values for the the endogenous variables
%   - params    [vector] vector of parameter values
%   - check     [scalar] set to 0 if steady state computation worked and to
%                    1 of not (allows to impose restrictions on parameters)

% read out parameters to access them with their name
% NumberOfParameters = M_.param_nbr;
% for ii = 1:NumberOfParameters
%   paramname = M_.param_names{ii};
%   eval([ paramname ' = M_.params(' int2str(ii) ');']);
% end
% 
% NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
% for ii = 1:NumberOfEndogenousVariables
%   varname = M_.endo_names{ii};
%   eval([varname '= ys(' int2str(ii) ');']);
% end

beta=NaN; %make parameter known to Matlab function, prevents crashes due to Matlab function with same name;
          %will be overwritten next
          
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = M_.param_names{ii};
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end

% read in instrument values
for ii = 1:size(options_.instruments,1)
  eval([options_.instruments{ii} ' = ys(strmatch(options_.instruments{ii},M_.endo_names,''exact'')) ;']);
end

% sa(1) = ys(strmatch('sa1',M_.exo_names,'exact'))
% sa(2) = ys(strmatch('sa2',M_.exo_names,'exact'))

% initialize indicator
check = 0;

%%%%%%%%%% Don't change the above code %%%%%%%%%%%%
% This code is for steady state of Ramsey with 1 sector, the gov tool is
% taux, we need to provide a steady state taking as given any taux, so
% treat taux as an exogenous/known variable (or a constant) here

%global theta dbar  L gL alpha rbest rho sigma ns n betas d alphas fspace_smolyak ilinear ilog_rL delta beta


global L n alpha  d


% private eqm
r_private = rbest;
n  = 2;
ns = 1;
ilog_rL = 1; % use exp(rL)/(1+exp(rL))

dij = dbar;
d=dij*ones(n); %dij is (1+iceberg) cost producing in j selling to i
d=d-(dij-1)*eye(n);

for i=1:n
    eval(['L(',num2str(i),',1)=L',num2str(i),';'])

    eval(['alpha(',num2str(i),',1)=alpha',num2str(i),';'])
end

% for i=1:n
%     alpha(i)=alpha(i)*exp(sa(i));
% end

% alpha1= alpha;
% alpha(1) = alpha1(1)*exp(0.01);
betas = 1;

%% solve for the steady state of gov assuming dG2/dT=0!!!
% This is not the true steady state, it's only for us to get initial guess
% --- steady state ---


options = optimoptions('fsolve', 'TolFun', 1e-12);
options.Display = 'off';

xguess(1:2) = 1.1;


[xval,fval_sol,eflag_ss] = fsolve(@(xval) solve_1s_2c_private_ss_nogL(xval,M_,ys),xguess,options); %iter

% guess wage and share of labor in each sector

w(1:n,1) = xval(1:n); % row vector
taf = 0;

rL(1,1) = exp(r1)*L(1);%rbest*L(1)
rL(2,1) = rbest*L(2);

T = alpha.*rL/delta1;

sextax = zeros(n);
for i=1:n-1
    sextax(i+1) = taux;
end
staf = zeros(n);

for i=1:n-1
    staf(i+1) = taf;
end

pi = zeros(n,n); % (1,2,3): country 2 export to ecountry 1 in sector 3, NOTE country order is opposite to our note
%d, matrix of n times n, dij country j ship to country i

xnimatrix1= T'.*((w'.*(1+sextax).*(1+staf).*d).^(-theta1));
 xn1 =sum(xnimatrix1,2).*ones(n,n);

pimatrix1 = xnimatrix1./xn1;
pi = pimatrix1;
xnimatrix= xnimatrix1;

betas =1;
x  =  (1+theta1)/theta1*w.*(L-sum(rL,2));

for m=2:n
    x(1)=x(1)+betas(1)*taux/(1+taux)*pi(m,1)*x(m);
end
tempt=0;

for m=2:n
    tempt=tempt+betas(1)*taf/(1+taf)*pi(1,m);
end

x(1)=x(1)/(1-tempt);

P=(sum(xnimatrix1,2)).^(-betas./theta1);

fval(1) = P(1)-1;
fval(2) = pi(1,2)*x(1)-pi(2,1)*x(2);

Tss = T;
xss = x;
Pss =P;
pi_ss=pi;


%(1+theta1*pi(2,2))/(1+taux)-theta1*pi(2,2)


xval_ss= xval;
wss = xval(1:n);
taf_ss = 0;
rss = rbest;


% construct M
% for i=1:n
%     GTss(i)= log(xss(i)^(-sigma)*Pss(i)^(sigma-1)*wss(i)/alpha(i));
%     Mss(i) = wss(i)/alpha(i);
% end

uc = xss(1)^(-sigmaH);
    
r2 = log(rbest);
%r1 = log(rbest);
rmul = 0;

%% write to our variable names in dynare
for i=1:n
    eval(['x',num2str(i),'=log(xss(',num2str(i),'));'])
    eval(['w',num2str(i),'=log(wss(',num2str(i),'));'])
    eval(['P',num2str(i),'=log(Pss(',num2str(i),'));'])
    eval(['riskfree',num2str(i),'=1/betad;']);

    eval(['T',num2str(i),'=log(Tss(',num2str(i),'));'])
    %eval(['r',num2str(i),'=log(rbest);'])
    eval(['sa',num2str(i),'=0;'])

    for ii=1:n
        eval(['pi',num2str(i),num2str(ii),'= pi_ss(',num2str(i),',',num2str(ii),');'])
    end
end

%% end own model equations, below from dynare

params=NaN(NumberOfParameters,1);
for iter = 1:length(M_.params) %update parameters set in the file
  eval([ 'params(' num2str(iter) ') = ' M_.param_names{iter} ';' ])
end

% NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
% for ii = 1:NumberOfEndogenousVariables
%   varname = M_.endo_names{ii};
%   eval([varname '=ys(' int2str(ii) ');']);
% end

NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
for ii = 1:NumberOfEndogenousVariables
  varname = M_.endo_names{ii};
  eval(['ys(' int2str(ii) ') = ' varname ';']);
end
